#Crevasses on moons 
#calculates fracture width
#robert.law@uib.no 2024

import os
import sys
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.integrate import quad

os.chdir(os.path.dirname(sys.argv[0]))
dir = os.path.dirname(os.path.realpath(__file__))

#equations from https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008GL036765 (auxillary information)

#variables
rhoi = 917 #(kg m-3) density of ice
rhow = 1000 #(kg m-3) density of water
g = 9.81 #(m s-2) gravity 
mu1 = 1.5e9 #(Pa) shear modulus 
mu2 = 0.5e9 #(Pa) second mu value
nu = 0.3 #Poissons ratio
Rxx = 0.1*1e6 #(Pa) background tensile stress
d = np.arange(50, 32500, 15) #(m) crevasse depth
#d = 2 * 1000 #(m) depth of crevasse
a = 0

nice_fig = True

#initialise output
int_out_mu1 = np.zeros_like(d)
int_out_mu2 = np.zeros_like(d)

#define G 
def G(b, d, g, mu, nu, Rxx, a):  
    
    alpha = 1 - nu
    omega = (d*d-b*b)**(1/2)
    psi = ((d*d-a*a))**(1/2)
    sigma = Rxx - ( (2*rhoi*g*d)/np.pi ) - rhow*g*a + ( (2*rhow*g*a)/np.pi )*np.arcsin(a/d) + ( (2*rhow*g*psi)/np.pi )

    A = (4*alpha*sigma/mu)*omega
    B = ( (4*alpha*rhoi*g*d)/(mu*np.pi) )*omega 
    C = ( (4*alpha*rhow*g)/(mu*np.pi) )*omega*psi
    D = ( (4*alpha*rhoi*g*b*b)/(2*mu*np.pi) )*np.log( np.abs((d + omega)/(d - omega)) )
    E = ( (4*alpha*rhow*g*(b**2 - a**2))/(2*mu*np.pi) )*np.log( np.abs((psi + omega)/(psi - omega)) )
    F = ( (4*alpha*rhow*g*b*a)/(mu*np.pi) )*np.log( np.abs((a*omega + b*psi)/(a*omega - b*psi)) )
    G = ( (4*alpha*rhow*g*a*a)/(2*mu*np.pi) )*np.log( np.abs((omega + psi)/(omega - psi)) )

    return A + B - C - D + E - F + G

mu = mu1
print(G(1, 22000, g, mu, nu, Rxx, a))
#sys.exit()

#for loop
for i in range(len(d)):
    #fix functions to values required
    print('Depth is: ' + str(d[i]) + ' m')

    int_fixed = lambda b: G(b, d[i], g, mu, nu, Rxx, a)

    #calculate integrals within bounds
    int, err = quad(int_fixed, 0, d[i])

    int_out_mu1[i] = int

mu = mu2 #new mu value
for i in range(len(d)):
    #fix functions to values required
    print('Depth is: ' + str(d[i]) + ' m')

    int_fixed = lambda b: G(b, d[i], g, mu, nu, Rxx, a)

    #calculate integrals within bounds
    int, err = quad(int_fixed, 0, d[i])

    int_out_mu2[i] = int

if nice_fig:
  
  fig, ax1 = plt.subplots(figsize=(80/25.4, 80/25.4))

  mpl.rcParams['font.size'] = 8

  r = ( ((int_out_mu1)/(np.pi)) )**(1/2) #calculate radius of sphere required, assuming fracture width matches 2/3 of sphere width
  #V1 = (4/3)*np.pi*r*r*r #volume 1
  V1 = ((4/3)*r)*int_out_mu1

  ax1.plot(d/1000, r/1000, color='black', lw=2)

  r = ( ((int_out_mu2)/(np.pi)) )**(1/2) #calculate radius of sphere required, assuming fracture width matches 2/3 of sphere width
  #V2 = (4/3)*np.pi*r*r*r #volume 2
  V2 = ((4/3)*r)*int_out_mu2

  ax1.plot(d/1000, r/1000, color='black', lw=2, linestyle='--')
  ax1.set_ylabel('Sphere radius ($km$)', fontsize = 8)
  ax1.set_xlabel('Fracture depth potential ($km$)', fontsize = 8)
  ax1.set_ylim(np.min(r)/1000, (np.max(r)/1000)*1.05)
  ax1.tick_params(axis='y', labelsize=8)
  ax1.tick_params(axis='x', labelsize=8)
  ax1.set_xlim(0,  32)

  ax2 = ax1.twinx()
  ax2.plot(d/1000, V1*1e-9, color=(0.76, 0.83, 0.89, 1), lw=2)
  ax2.plot(d/1000, V2*1e-9, color=(0.76, 0.83, 0.89, 1), lw=2, linestyle='--')
  ax2.set_ylim(np.min(V1)*1e-9, (np.max(V2)*1e-9)*1.1)
  ax2.set_ylabel('Sphere volume ($km^{3}$)', color=(92/255, 142/255, 181/255, 1), fontsize = 8)
  ax2.tick_params(axis='y', labelsize=8)
  

  plt.tight_layout()
  plt.savefig('Sphere_volume.png', dpi = 400)

  plt.show()
    
else:

  #plot Fig. 3
  #plt.plot((400*int_out)**(1/2), d) #this one for lake diameter, not sphere diameter
  r = ( ((6*int_out)/(4*np.pi)) )**(1/2)

  plt.plot(d, r/1000, color='black', lw=2)
  plt.ylabel('Sphere radius ($km$)')
  plt.xlabel('Fracture depth potential ($km$)')
  plt.ylim(np.min(r)/1000, (np.max(r)/1000)*1.05)
  plt.xlim(0,  np.max(d))

  #plt.axhline(y=980, color='black', linestyle='--')

  plt.show()

sys.exit()

#plot Fig. 1
plt.plot(d, int_out)

plt.xlabel('d')
plt.ylabel('fracture area')
plt.xlim(0, 250)
plt.ylim(-4, 12)
plt.show()





sys.exit()
#old (non-looped)

#fix functions to values required
int_fixed = lambda b: G(b, d, g, mu, nu, Rxx, a)

#calculate integrals within bounds
int, err = quad(int_fixed, 0, d)

print(int*1e-9)




